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Abstract 



We study systematically stationary solutions to the coupled Vlasov and Poisson equa- 
tions which have 'self-similar' or scaling symmetry in phase space. In particular, we 
find analytically all spherically symmetric distribution functions where the mass density 
and gravitational potential are strict power laws in r, the distance from the symmetry 
point. We treat as special cases, systems built from purely radial orbits and systems 
that are isotropic in velocity space. We then discuss systems with arbitrary velocity 
space anisotropy finding a new and very general class of distribution functions. These 
distributions may prove useful in modelling galaxies. Distribution functions in cylin- 
drical and planar geometries are also discussed. Finally, we study spatially spheroidal 
systems that again exhibit strict power-law behaviour for the density and potential and 
find results in agreement with results published recently. 
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1 Introduction 



Star clusters, dark matter galactic halos, and clusters of galaxies are essentially col- 
lisionless self-gravitating systems obeying the coupled Poisson-Vlasov equations and 
therefore equilibrium solutions to these equations are of great importance. While sub- 
stantial progress has been made through numerical simulation, there has also been a 
persistent school of analysis that seeks to obtain exact, analytic results, particularly in 
asymptotic or 'stationary' limits. 

One branch of this latter school has studied the evolution of coUisionless matter in an 
expanding universe (Fillmore & Goldreich 1984; Bertschinger 1985; Gurevich Sz Zybin 
1988, 1990; Ryden 1993). The solutions are time dependent of necessity but show steady 
or at least adiabatic behaviour at late times. The treatments make use of an intuitive 
self-similar symmetry which can seem rather ad hoc (even if exceedingly clever) and so, 
difficult to access and generalize. Moreover, there are some remaining inconsistencies 
in the various published results. Our intention here is to study such coUisionless self- 
similarity in a simple and systematic way, beginning with strictly stationary examples. 

Self-similar symmetry has also been found in the study of coUisional systems such as 
the cores of globular clusters (e.g. Lynden-Bell 1967; Lynden-Bell and Eggleton 1980; 
Inagaki and Lynden-Bell 1983, 1990). These systems allow the study of the evolution 
towards core-halo configurations but never yield true thermodynamic equilibria in the 
form say of stationary power law behaviour. Nevertheless in various intermediate (spa- 
tially and temporally) stages simple power law do appear and may in fact correspond to 
the solutions we find for stationary coUisionless systems. It is well known in the hydro- 
dynamic literature (e.g. Barenblatt and Zel'dovitch 1972, hereafter BZ) that self-similar 
solutions arise as intermediate asymptotics between boundaries, and it is probably for 
this reason that they are found in the pre and post collapse phases of the coUisional 
systems. One might speculate that they are as close to 'equilibrium' as such systems 
get. In addition to globular clusters, these systems may arise in the intermediate stages 
of collapse of a star cluster to a black hole and in the ultimate state of coUisionless dark 
matter halos. 

In this paper we study systematically the family of distributions that are exactly station- 
ary and that possess both precise geometric (usually spherical) and scaling symmetries. 
These distributions are solutions to the coupled Vlasov and Poisson equations: 



Here / is the usual mass distribution function, $ is the gravitational potential, H 
is the Hamiltonian and canonical pairs. The solutions we find include the 

intermediate asymptotic limits of the relevant preceding studies. They overlap directly 




(1.1) 




(1.2) 
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with the recent work of Evans (1994, hereafter E94) on "power law galaxies" and to a 
lesser extent with the models of Tremaine et al. (1994). These latter models are 
similar in spirit to our own but are not strictly scale free and so prove to be rather more 
complicated to express and use. Nevertheless they share the observationally important 
property of a cusped central density and possess the substantial advantage of having 
finite mass. 

The power-law galaxies of E94 (see also Evans and de Zeeuw 1994, hereafter EdZ) are 
scale invariant and therefore simple to use. They were constructed to provide a family 
of versatile models of galaxies that could be used to analyse observable properties of 
galaxies such as line profiles from stellar absorption-line spectra. Our approach is more 
systematic and in fact, for the case of spherical symmetry with velocity space anisotropy, 
yields a more general class of solutions. We mostly leave applications such as those found 
in van der Marel Sz Franx (1993), E94 and EdZ for future publications. 

We also use this paper to introduce to the astronomical community a systematic method 
of classifying scaling or self-similar symmetry first given by Carter and Henriksen (1991, 
hereafter CH) in a rather formal style. The technique is easy to use and yields the 
reduced equations (with manifest scaling symmetry) in a standard form. In addition, 
the method includes the various scaling symmetries that produce irrational power laws 
formerly referred to as scaling symmetries of the 'second' kind (e.g. BZ). 



2 Radial Orbits in Spherical Symmetry; 

An Introduction to the CH Analysis of the Vlasov-Poisson Equations 

We introduce the canonical distribution function F where 

fM^%^Hvo)SM . (2.1) 
The Vlasov and Poisson equations become 

VrdrF - dr^dy^F = , (2.2) 



drir^dr^) j Fdvr . (2.3) 

Following CH we seek a self-similar or scaling symmetry in phase space by requiring 
that the distribution function satisfy the equation 

UF = (2.4) 
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where 



>Ck = k^dj = 5rdr + Wrd^^ (2.5) 

is the Lie derivative with respect to the (phase space) vector operator k. In other words, 

k is the scahng direction in phase space. It is convenient in these formulae to imagine 
that r and are scaled interms of fiducial values (not to be confused with real constant 
lengths). 

Equation (2.4) holds so long as F = F{Q) where C, = r^^l^^ jvr- The dimensionless 
real number bjv gives the similarity 'class' of the symmetry in the sense of CH and 
it is generally fixed only by boundary conditions or by the dimensions of conserved 
quantities. 

Our first step is to choose new phase space coordinates to replace r and Vr- Following 
CH we define a new 'radial' coordinate R{t) such that 



£k = k'dj = Or . (2.6) 



k^djR = 1 and we obtain the transformation laws 



r I 6 



(2.7) 



and 



^ = e-^'^sgniS) . (2.8) 

The orthogonal invariant or self-similar variable X (in the notation of CH) satisfies 
k^djX — 0. In this simple example X is an arbitrary function of (. Without loss of 
generality, we take it to be 



X =\S l-^^/-^) /C = Vre-""^. (2.9) 

As of yet, there are no restrictions on the solutions; (X, R) is merely a new coordinate 
system. The restrictions appear only subsequently when we impose the invariance of 
various quantities under the action of k. 

While the scaling of dimensional quantities under the action k is uniquely determined, 
the form is hardly standard and varies from author to author. CH present an alternative 
analysis introducing a 'dimensional' algebra in an appropriate 'dimension space'. In the 
current example it is sufficient to choose a dimension space consisting of length, velocity 
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and mass. There is then a three parameter multiphcative rescahng group with elements 
A = (6*^,6^, e'^) that describes the scahngs of dimensional quantities. Alternatively, if 
we consider changes in the logarithms of these quantities, we have an additive rescaling 
group with elements a = (5, z^, //) . The vector components of A or a correspond to the 
scaling in length, velocity and mass respectively. 

Each dimensional quantity ^' in the problem has its dimensions represented by a di- 
mensionality (co)vector in the dimension space, and the change in the logarithm of 
the quantity is given by (CH) 

= Or-^ = (dvE, • a)* . (2.10) 

Strictly speaking, the vector k should correspond to the rescaling vector a. In other 
words, we should replace equation (2.5) with 

k-'dj = Srdr + vvrdy^ + jimdm . (2-11) 

However as we will soon see, the invariance of G under rescaling implies that mass 
rescaling is not independent of length and velocity rescaling. This allows us to reduce 
our scaling algebra element to a = (5, u). 

The dimensional quantities in the current problem F, $ and G have the following 
dimensionality covectors in the chosen dimension space (length, velocity, mass); 



d^ = (-1,-1,1) , 

cl$ = (0,2,0) , (2.12) 
dG = (1,2,-1). 

The requirement that G be invariant under the rescaling group action implies a -de = 0, 
or on performing the scalar product (direct multiplication since d is a covector) 

li = 2u + 5. (2.13) 

Consequently the dimension space may be reduced to the sub-space of {length, velocity) 
wherein the rescaling group element becomes a = (5, u) and 



dF = (0,l) , 

d$ = (0,2). (2.14) 
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The sub-space algebra element a now corresponds to our choice of the scaling vector k. 

The scaling symmetry can be imposed on F and $ by equation (2.10) which, with 
equation (2.14) requires 



F{X,R) = F{X)e''^ , 

$(X,i?) = ¥(X)e^^'^ . (2.15) 

These equations can be simplified further by noting that at fixed r the potential is 
independent of Vr or equivalently 

$ = constant . (2.16) 
Substituting into the Vlasov and Poisson equations yields 

dlnF _ 

dlnX ~ X^ + 2^ ' 

and 



(2.17) 



2y- + 2j^ = G\^-j J FdX . (2.18) 

It is now plain that the class of solutions depends only on the ratio 6/iy and so without 
loss of generality, we set = 1. 

Equation (2.17) now integrates to yield F — C \ X"^ + 2$ and therefore by 

equations (2.9) and (2.15) 

F = C\v^ + 2^ = C\2E . (2.19) 

Here C is a normalization constant given below. The potential is given by $ = $e^"^ = 
$ (I 5 I r)'^^^ . It is convenient to introduce an explicit fiducial length a. The potential 
can then be written as 



/ J. X 2/(5 

$ = (-) (2.20) 



where $o has units of velocity^ and replaces $. We also find from (2.18) that 
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For bound solutions, $ (or equivalently must be negative and therefore d < —2. 
Formally 5 > and $ > also give an inwardly directed gravitational force, but now 
there is no natural cut-oflF in equation (2.18) and the whole integral does not converge. 

The mass density obeys a simple power law: 



pocr(2/^-2), (2.22) 

and 6 < —2 implies that the density law lies in the range — 2 > | — 2 > —3. In addition, 
the velocity moments of the distribution function also follow simple power laws: 



y2n 0^ ^2n/5 ^2.23) 

indicating that velocities increase without limit as one moves towards the center of the 
distribution. 

It is interesting that this distribution function has the same functional dependence on 
energy as a Plummer model with n = 5/2. However, the density does not satisfy 
the Lane-Emden equation since it is 'cusped' at the centre and so does not have the 
right boundary conditions (e.g. Binney & Tremaine 1987; hereafter BT). It appears 
to be the distribution function corresponding to the asymptotic solutions found by 
Fillmore and Goldreich (1984). Indeed, although their work was done in the context 
of an expanding universe, their solutions in the case of spherical symmetry turn out 
to be time- independent (see their equation (39)). Figure 1 shows a contour map of 
the distribution function for S = —17. This value of S is chosen so that our Figure 1 
corresponds to their Figure 10. Contours give equal probability density in the correct 
phase space and so the integral of the contour value times the area element for any 
region of the plot gives the total number of particles in that region of phase space. 
In addition particle orbits are coincident with these contours and in fact Fillmore and 
Goldreich (1984) construct their plots by calculating particle orbits. 

We are left with the free parameter 5 in the range (— 2, — oo). How might this be 
determined? It is clear from the discussions of earlier papers that in general it appears 
as a sort of eigenvalue when the solution is seen to arise from more general initial 
conditions. However another manner in which it may be fixed is to require additional 
global invariants. For example, if we introduce a characteristic mass, which may be 
either the mass of each particle or that of a central massive object (or both in the sense 
that once there is one fixed mass we can measure all others in terms of it). Then there 
will be no mass scaling in the scaling algebra which means fi = 0. From this and equation 
(2.13) we see that 6 — —2 \ This corresponds to a central point mass surrounded by 
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massless particles that nevertheless are distributed in energy space according to equation 
(2.19). 

Having illustrated our method in detail in this simple example with radial orbits, we 
proceed to give our results briefly for other cases of interest. 



3 Isotropic Orbits in Spherical Symmetry 

Here, the full 3D distribution function has the functional form / = f{r,v), where 
v'^ = + Vq + v'^. The Vlasov and Poisson equations are now respectively {vr ^ 0) 



Vdrf - idr^)dj = 0, (3.1) 

^drir^dr^) = AttG J ^Wfdv. (3.2) 

The analysis proceeds as above and in fact equations (2.6), (2.7) and (2.8) apply here 
as well while equations (2.5) and (2.9) hold with replaced by v. The dimension 
space is again (length, velocity, mass). The dimensionality covector of / becomes d/ = 
(—3, —3, 1) while those of G and ^» remain as in equation (2.12). The invariance of G 
requires equation (2.13) as before so that the dimensionality covectors in the subspace 
of (length., velocity) are 



d/ = (-2,-1) , 

d$ = (0,2) , (3.3) 



whence following equation (2.10) 



/ = 7(^)e-(^^+^)^ , 

$ = $(X)e2^^ . (3.4) 

We will set 1/ = 1 in the sequel to be consistent with the notation in Section 2. Again, 
things simplify as $ = constant. Substituting into equations (3.1) and (3.2) yields 

dlnf ^ (2r) + l)X 

dX X2 + 2$ ^ ' ' 

and 
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_ ^a/-2<I> _ 

2(5 + 2)$ = (47r)2G / X^f{X)dX. (3.6) 

Equation (3.5) easily integrates and we find 

7 = C I + 2$ _ (3_7) 

It becomes clear from this last equation that for a bound system ($ < 0) once again 
5 < —2. In this case however, unbound solutions do exist since the integral in (3.6) 
converges for 5 > 1. 

We may use equations (3.4) and (2.9) (with v in place of Vr) to write the solution for 
the distribution function as 



/ = C\v^ + 2^ |-(m/2)= c I 2E |-('5+V2)^ (3.8) 
and that for the potential as 

^ = ^{\6\rf^ = ^a{iy^'. (3.9) 

As before, a is a fiducial length and $a is a constant with units velocity'^ that replaces 
The normalization constant C is determined from the Poisson equation: 

_ f -{l + S/2){l-6)T{l-6) \ {-2^a)' 

\ 27rV252r (1/2 - 5) / Ga2 " ^'^'^^^ 

It is interesting to note that the corresponding density law has the same functional form 
as in the pure radial case discussed above. Moreover, f ^"^ oc r^"^/*^. 

Once again, at the cost of admitting a singular potential and density we have found a 
large class of isotropic analytic solutions. The distribution function of equations (3.8) 
and (3.9) is much like the Plummer model (e.g. BT, p223), but the potential member 
of the pair is singular. The stability theorems of Antonov and of Doremus, Feix and 
Baumann (as cited in BT section 5.2) seem to imply moreover that this distribution 
is stable. The velocity dispersion for both radial and isotropic cases increases towards 
towards the centre of the distribution even though the mass at the centre is zero. (The 
mass inside radius r diverges r — > cxo in common with the standard 'isothermal' model.) 
These solutions therefore might be useful models for star clusters at the centers of active 
galaxies. For example, if a detailed mapping of the velocity dispersion in M87 were to 
show a slower than Keplerian increase with decreasing r, then star cluster models of 
this type might provide an alternative to the massive black hole theory. 
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We note that once again the hmit 5 = —2 corresponds to an invariant central mass 
surrounded by massless particles. This is the Keplerian, and possibly the central black 
hole, limit. 



4 Systems with Cylindrical and Planar Symmetry 

It is straightfoward to extend the previous analysis to systems with cylindrical and 
planar symmetry. Such solutions might provide insight into the filaments and sheets 
seen both in N-body simulations and in large-scale galaxy surveys. Our solutions should 
correspond to those found in Fillmore and Goldreich (1984) though our solutions are 
stationary and do not include cosmological expansion. 

For systems with cylindrical symmetry, we choose the canonical distribution function 

/(r,v) = ^^^^^^Kve)Kv.) (4.1) 

where w is the distance from the axis of symmetry. The Vlasov and Poisson equations 
are then 



d^{wd^^) = G j Fdv^ . (4.3) 
We find the family of solutions 

/ oc $ oc w""/' . (4.4) 

The mass /length inside a radius w is iJi{w) oc w'^^^ while the density is oc ■uj'^l^~'^ . If we 
require that the density decrease with increasing w (but allow for a cusped density law 
at the origin) and also require that there be no central mass concentration (//(tu) 
for tt7 — > 0) then 8 is constrained to be greater than 1. 

For planar symmetry, we choose 

j{v,^) = F{z,v,)b{v,)b{vy) . (4.5) 
The Vlasov and Poisson equations are then 

v,d^F-d,^d,F^^ , (4.6) 
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(4.7) 



In this case, the family of solutions is 



/ oc $ oc z^" . (4.8) 

1 < 5 < 2 is enough to insure that the density will decrease with increasing z and the 
mass per unity area will vanish for 2; — > 0. 



5 Anisotropic Orbits in Spherical Symmetry 



This example is more challenging then the others treated in this paper and the re- 
sults are somewhat surprising. The distribution function depends on three phase space 
coordinates {r,Vr,j'^) where p = r^(f^ + f|) is the square of the transverse angular 
momentum (Fujiwara 1983). The Vlasov and Poisson equations become respectively 



Vrdrf + ( ^ - 5,$ ) {d,J) = 0, (5.1) 

drir^dr^) -Att'^G I df I f{r,Vr,f)dvr = 0. (5.2) 



The dimension space is taken to be (length, velocity , (angular momentum)'^ , mass), with 
the scaling algebra vector a = {S,u, X,fi). However the invariance of G allows us to 
work in the reduced scaling space (length, velocity, (angular momentump ) . Moreover 
we know that we can set one of the scale factors equal to unity since only the ratios 
have physical significance, and we choose this to be = 1 as above. The symmetry we 
seek can therefore be written in the simplified form 



¥dj = rSdr + Vrdy^ + Xfdp . (5.3) 

As before we choose R to lie along this direction so that equations (2.6), (2.7) and (2.8) 
continue to apply. 

There must now be two invariant coordinates X* orthogonal to R which satisfy k^djX'^ = 
0. These are linear, first-order partial differential equations that integrate easily to give 

X' =T{(:) xvre-^, (5.4) 



11 



where are arbitrary functions oi C, = Vj. j It is convenient to choose X^^^ = 1 
and X^^^ = 1/C so that 

X«=X = t;^e-^ (5.5) 

and 

We shall write our equations in terms of R, X, and Y. 

In the reduced dimension space the quantities have the dimensionality covectors 

d/ = (0,1,-1), 

d$ = (0,2,0), (5.7) 

whence as usual 

/ = /(x,y)e(i-^)^, 

$ = $(X,y)e(2^^ . (5.8) 

The potential at fixed r should be independent of Vr and and we impose $ = constant. 
The Vlasov and Poisson equations now become 

-{26 + 1)7 - Xdxl - 2(1 + S)YdyJ + -^{S^Y - 2^)Xdxf = 0, (5.9) 
2$ Q + ^) - 47r^G' J dV J dXj{X, Y) = 0. (5.10) 

where 



A = 2(1 + 5) (5.11) 

follows from the requirement that these equations be independent of R. We see from 
equation (5.10) that S < —2 for the system to be bound. 

Equation (5.9) is a quasi-linear, first-order partial differential equation that can be 
solved exactly. The characteristic equations are 

df dY XdX 

•' (5.12) 



(25 + 1)/ -2(1 + 6)Y {S^Y-X^-2^) ' 
the first of which is clearly integrable and gives 

_ { 2S + 1 \ 

f = F{^)Y \^^), (5.13) 
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where ^ is a function constant on the integral curves of the second equation of (5.12). 
This second equation may be integrated by introducing the variables u, s such that 



s = ^2 + 2$ , (5.14) 

and 



5^Y = u + s . (5.15) 

The equation to be integrated is now reduced to 

du 

u— + (2 + S)u+ (l + 6)s = , (5.16) 
as 

which is of a type already known to Leibnitz in 1691. The solution is immediate for 
s 7^ by changing the dependent variable to say y = u/s. We find thereby that the 
quantity ^ {X, Y) may be taken to be 

^ = y-'^ \ 5^Y + X'^ + 2^ (5.17) 

The solution is now determined but for the choice of the arbitrary function F{^). The 
constant $ is determined from equation (5.10) or equivalently 

_/-! 2\ /•-2¥/5^ i-\/-2W-S'^Y _ 

2^ T + T9 H^tt'C / dY ^ dXf. (5.18) 

In this formula / is used in the form of equation (5.13) and the order of integration may 
be reversed if the integrals exist. 

It is instructive to use equations (5.5), (5.6) and (5.8) to return to physical variables 
whence 



^ = j-2 I 2$ + -i;2 |'5+i= j-2 \2E\^+\ (5.19) 

and 

/ = j-(^)F(0 (5.20) 

where v'^ = + j'^/r'^ is the squared magnitude of the velocity. It is remarkable that 
despite the arbitrary nature of F(^), equations (2.20) and (2.22) continue to give the 
potential and density power laws. The amplitude $ does depend however on F{^) 
through equation (5.18) and one must insure that the integrals in this equation exist. 

Equation (5.20) represents a new class of distribution functions with velocity space 
anisotropy. Models with F{^) oc correspond to the power-law spheres discussed in 
E94. For these models, the distribution function takes the simple form 
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/ = Cf I E 



(5.21) 



where a and b are related to a and 5 in a simple way. (Actually, distributions with energy 
and angular momentum dependence given by equation (5.21) were first discussed by 
Camm (1952). There, the spatial dependence of the mass density and potential are given 
by solutions to the Emden-Fowler equation with appropriate boundary conditions.) As 
discussed by E94, it is straightfoward to calculate velocity moments for these distribution 
fvmctions^ For example, one can calculate the velocity space anisotropy parameter 



Here, P = {a = — {26 + 1) /2 (5 + 1)) corresponds to an isotropic distribution in ve- 
locity space, P = 1 (a =1/2 (5 + 1)) corresponds to purely radial orbits, and (3 ^ —oo 
[a — oo) corresponds to a distribution of purely circular orbits. Evidently, equation 
(5.21) describes a two-parameter family of distribution functions where one parameter 
specifies the density law and the other parameter specifies the anisotropy in velocity 
space. The distribution given in equation (5.20) is more general than this indicating 
that entirely diflFerent distribution functions can have the same density law and velocity 
space anisotropy. We illustrate this through the sequence of figures 2(a)-2(d). These 
figures show contour plots of distribution functions in velocity space ((j/r, Vr) space) 
for fixed r. Figs 2(a)-2(c) are power-law spheres described by equation (5.21) and corre- 
spond to an isotropic distribution (/? = 0); a distribution constructed from nearly radial 
orbits (/3 = 0.9); and a distribution constructed from nearly circular orbits (/? = —9.0) 
respectively. Fig 2(d) is a composite model constructed from the power-law models used 
in Fig 2(b) and 2(c), i.e., from nearly radial and nearly circular oribits. By construc- 
tion, (3 — as in the isotropic case though the actual distribution function is entirely 
different. 

It is interesting to note that the degeneracy of models discussed here (density law fixed 
while the functional form (admittedly of a specified argument in j and E) free) is in the 
opposite sense to the degeneracy found in the 'Einstein model' where particles move on 
all possible circular orbits in spherical symmetry while the density and potential laws 
are arbritrary (Fridman and Polyachenko 1984). 

It is plain that a great deal remains to be discovered about these distribution functions, 
which have the character of a 'post core collapse' (zero central mass) stationary state. 
For example, for which choices of F{^) are they stable? Fridman and Polyachenko (1984) 
show that the Camm models are linearly unstable for sufficiently large anisotropy, but 
are stable towards the isotropic limit. But linear stability is not the same as non-linear 
stability which has really to do with the existence or otherwise of asymptotic distribu- 
tions towards which / tends. By studying the stability of our stationary solutions we 
might discover a hint as to how to define a useful 'entropy' function that characterizes 



P = 



26 + 1 



(5.22) 



2(5 + 1) 



+ a 
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such asymptotic equilibria. This latter question has recently been given new life by the 
studies of TVcmaine et al. (1986), Wiechen et al. (1988) and Aly (1989, 1993). This last 
author has shown for example that the 'softened' Plummer model with 6 = —4 in our 
notation actually attains the minimum energy subject to a fixed mass and a fixed value 
of the 'entropy'. If this value is used naively in our models, one obtains p oc r~^/^ and 
$ oc r~^/^. Finally one might ask how closely can a stable model approach Maxwellian 
type distributions. 



6 Axially Symmetric Solutions with Ellipsoidal 
and Hyperboloidal Symmetries 



In order to make contact with the recent studies of E94 and of EdZ we show here how 

our method yields the strictly scaling subset of their axially symmetric solutions in a 
very direct way. The scaling symmetry in phase space for this case is actually simpler 
than that for the anisotropic spherical geometry dealt with above. We use cylindrical 
coordinates w,(f), z and write the Vlasov equation in the symmetric form 

Vr.d^ f + v,d, f+('^-d^^]d,^f- ^a,, / - / = 0, (6.1) 

while the Poisson equation becomes 

-d^{wd^^) + {d^f^ = AttG I fdv^dv^dv, . (6.2) 
w J 

Following E94 we seek solutions with confocal ellipsoidal symmetry wherein $ = $(tt) 
where 

^2 

u^ = w^ + ^ . (6.3) 

We differ from E94 in that we do not introduce a 'core' radius into the problem, since 
strictly speaking this would prohibit the existence of true geometrically scaling solutions. 
It will be clear however from our method that one can add a core radius squared to 
without changing the solution, provided that it is ignored in the scaling symmetry. 
Thus our solutions as well as those of E94 do not 'really' know of its existence and so 
we prefer to suppress it. It is amusing to note also that the substitutions — > —q'^ and 
u ^ V where 

q^ 

in our results will convert them to hyperboloidal symmetry in which the equipotential 
surfaces are confocal hyperboloids of one sheet rotated about the z axis. Since these 
extend to infinity it may be that they are of little interest, but we note that near z — 
they tend to a configuration that might be associated with a vertically stratified disc 
containing a central 'hole'. 
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We proceed with a direct substitution of $ = into the right hand side of (6.2) 

followed by a collection of terms in equal powers of w to discover the necessity for a 
distribution function in the ansatz form 

/ = fi{u,v^,v^,Vz) + f2{u,v.^,v^,Vz)w'^. (6.5) 

Hence the Poisson equation splits into the two equations 

f 2$' 

47rG / fi dv^ dvA, dvz = 1 — tt, (6.6a) 

J u 

AttG J f2dv^dv^dvz = - 1^ - , (6.66) 
where the primes denote total derivatives with respect to u. 

The key procedure is the substitution of the ansatz (6.5) into the Vlasov equation 
(6.1) and the subsequent need to satisfy the equation for the coefficients of each of the 
powers of w and of z. Each such coefficient presents us with simple quasi-linear partial 
differential equations to solve and despite the apparent danger of over determining the 
problem, everything in fact works smoothly. Thus from setting the coefficient of z equal 
to zero we learn that 

/i = fi + ^''^^''^'^ 
/2 = /2(^ + $,^^,^<^) . (6.7) 
Then by setting the coefficient of equal to zero there foUows finaUy for /i 

v'^ -L t;^ -|- t;^ 

/i = /i I + ^ I ^ ME). (6.8) 

The vanishing of the coefficient of w yields then directly for /2 



h = viF^iE), (6.9) 

after which the coefficient of the last term in vanishes identically. Consequently we 
concur with E94 that the necessary ansatz (6.5) is in fact of the form 

/ = h{E)+w^vlF2{E) ^ h{E)+jlF2{E) , (6.10) 

although at present /i and F2 are arbitrary functions of the energy E. The Vlasov 
equation is now identically satisfied, in accordance with Jeans' theorem for stationary 
solutions. 



16 



As in the previous sections we turn now to the imposition of a scahng symmetry in 
phase space. The Lie derivative will be along the vector 



k^dj = u5du-\-v^dy^-\-v^dy^-\-Vzdy^ = da, (6-11) 
where R = R{u) and so as usual 

u=—sgn{5), (6.12) 

^ = e-^^sgn{5). (6.13) 

Proceeding exactly as before to solve for the invariants X^'^^ from k^djX^'^^ = gives 
the convenient choices 



(6.14) 



Moreover starting once again with the dimension space covector a = [5, u, //), preserving 
G (but not any characteristic length) and finding the dimension vectors of /i, F2 and 
$ in the reduced dimension space gives 

h = A(x(^))e-(^^+^)^, 

F2 = F2(xW)e-(^'^+3^)^, (6.15) 
$ = ¥(xW)e2'^-R. 



We find subsequently that u = 1 in accordance with the form assumed above for k. To 
proceed we first observe that by the scalings (6.14) and (6.15) 

E = ( ^i!^ + $ ] = :E(X(^))e2^. (6.16) 



Then the compatibility of equations (6.8) and (6.15) (together with the isotropic char- 
acter of /i in velocity space) requires 

h{E)=J,{Ee-'^)e-^''+'^^ . (6.17) 

It is clear therefore that fi{E) should be an homogeneous function of order a (i.e. 
fi{kx) — k"fi{x)) where 
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a = -(5+1/2). 



(6.18) 



In one dimension such a homogeneous function is a power law of power a so that 
f xix) — Ax'^ and consequently 



/i(^) = A£;«. (6.19) 

An exactly similar argument based on the compatibility of equations (6.9) and (6.15) 
shows that 

F2{E) = BEf^, (6.20) 

where 

/5 = -(25 + 3/2). (6.21) 

In the preceeding formulae, A and B are arbitrary constants and E should be interpreted 
as the modulus of E when E < 0. 

Thus imposing a strict scaling symmetry (which will of course be communicated to the 
particle orbits through the Hamiltonian by virtue of the scaling in equation (6.16)) leads 
in a self-contained fashion to the distribution function 

/ = 7i(x(^)) e-(2^+i)« + f,F2{X^^^) e-("^+3)^, (6.22) 
= AE'^ + Bj^E^. 

When proper account is taken of the scaling law of specific angular momentum (oc 
Q-iS+i)R |-,y usual dimension space arguments) equation (6.21) shows explicitly that 
both terms in equation (6.22) scale in the same way, as they should in order that / have 
this same strict symmetry. 

We note that the potential varies as $ oc ti^/*^, which on comparing with equation (2.1) of 
E94 shows that his parameter (3e related to our scaling parameter 6 by Pe = ~(2/5). 
Our result (6.22) is then equivalent to the form (2.6) of EdZ when the core radius is 
zero, as it should strictly be. 

All of the discussion of this case is now reduced to the quadratures in equations (6.6) 
which, together with (6.22) and the scaling symmetry of equations (6.15), become 

AttGA J E^'dX^'Ux^'Ux^'^ = $(^45+^^^^) , (6.23) 
AttGB J E^X^^^fdX^^Ux^^Ux^^^ = 4¥ (^-l) ^'(<^-l). (6.24) 

Prom equation (6.23) one observes that for bound solutions $ < 0, we require 

5 < sgn{A) , (6.25) 
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and by equation (6.24) 



\ -1] (d -l)sgn{B) < 0. (6.26) 



If both A and B are presumed positive then there is no cutoff scale except the zero 
energy surface. However there seems to be no need for this to be true, and in general 
one can have an energy cutoff for any and vice versa. For unbound solutions with E 
and $ both positive, we note that 6 should be positive if the gravitational acceleration 
is to be directed inwards. 

The evaluation of the integrals (6.24) and (6.25) has been extensively discussed by E94 

and EdZ and are easily done by changing variables to E,6 = tan~^ I f + j /vg, 

and (f) = tan~^ {V(i)/v^). In addition, one can calculate various moments of the distribu- 
tion function and for example, derive theoretical line profiles for stellar absorption line 
spectra. 



7 Conclusions 



In this article we have introduced a systematic method of imposing scaling symme- 
tries or self-similarity in phase space. These symmetries, along with various geometric 
symmetries, allow one to find a wide range of analytic solutions to the coupled Vlasov 
and Poisson equations. This has led us to recover in a straightforward way the known 
'power law' or scaling solutions, some of which are just beginning to be recognized for 
their practical implications. In addition, we have found a new and very general class 
of spherical systems with velocity-space anisotroy. The method is sufficiently powerful 
that we expect it to be useful in time dependent cases as well. At the very least we will 
be able to reduce these problems by one variable as was achieved here for example in 
section 5. 

The self-similar solutions discussed in this paper have applications in a wide variety of 
astrophysical problems. For example, the solutions found in sections 3 and 4 provide 
examples of distribution functions where the velocity dispersion increases without limit 
as r — (as per equation 2.23). These models might be relevant to observations of 
galactic nuclei where the densities are apparently correspondingly 'cusped' (Tremaine 
et al. 1994). Although this might be thought to cast doubt on various black hole 
'detections', it must be remembered that the collisionless constraint has less and less 
relevance on galactic time scales in these dense regions. Thus the system will certainly 
evolve and be non-linearly unstable under the influence of these collisions. Nevertheless 
it is conceivable that recently formed structures might show this behaviour. One ap- 
plication where the system is almost certainly collisionless, even over a Hubble time, is 
the collapse of the dark matter. Here, collapse can occur along one, two, or three axis 
and so, following Fillmore and Goldreich (1984) we have studied self-similar collapse in 
planar and cylindrical geometries. 
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In section 5 we have found the power-law solutions that are spherically symmetric though 
anisotropic in velocity space. A subset of these models are the 'power-law galaxies' of 
E94. However, our models are more general as illustrated in Figures 2(a)-2(d). The 
models therefore allow for more freedom in modelling galaxies and fitting observables 
such as absorption-line spectra to theoretical predictions. 

Finally section 6 has enabled us to demonstrate the ease with which the spheroidal so- 
lutions may be found when the spatial symmetry and the scaling symmetry are imposed 
separately. Our conclusions are the same as those of E94 and of EdZ94, although we 
do point out that the core radius is not actually playing a role in these solutions and 
it is apparent that these solutions form part of the wider family of solutions discussed 
in this paper. We also remark that hyperboloidal solutions of this type obviously also 
exist. 

Future work will focus on the time dependent solutions with and without an expanding 
background partly in hopes of obtaining the results of Fillmore and Goldreich (ibid) 
without the singularities in the density, but also in order to make a survey similar 
to the present survey of the stationary solutions. More detailed investigation of the 
stability and other properties of the solution in section 5 is also required. 
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Figure Captions 



Figure 1. Contour plot of the distribution function for a spherically symmetric system 
with radial orbits. The distribution function is given by equations (2.19) and (2.20) with 
5 = —17 (chosen so that this figure corresponds to Fig. 10 of Fillmore and Goldreich 
(1984)). Horizontal and vertical axis are r (measured in units of some fiducial length 

1/2 

'a' and (measured in units of ). Contours are linearly spaced and range from 
0.04 to 1.0. 

Figure 2(a). Contour plot of the distribution function for a system with that is 
isotropic in velocity space (/? = 0). The distribution function is given by equation 
(5.21) or equivalently (5.20) with F(^) = For definiteness, we take 5 = —8 and 

a = —1.07. r is fixed and the horizontal and vertical axis are [vq + v"^) = j/r and 

Vr measured in units of ^>(r)^/^. The contours are logarithmic ranging from 1 to 10^ 
with the distribution function increasing towards the lower left of the plot. (We have 
not worried about normalization of the distribution function.) 

Figure 2(b). Same as Fig. 2(a) but for a distribution function constructed from nearly 
radial orbits (/3 = 0.9, S = -8., a = -0.17). 

Figure 2(c). Same as Fig. 2(a) but for a distribution function constructed from nearly 
circular orbits (/? = -9.0, S = -8., a = -10.07). 

Figure 2(d). Same as Fig. 2(d) but here the distribution function is given by the sum 
of the two power-law distributions used in Figs 2(b) and 2(c). By construction, this 
distribution has (3 — 0. 
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